phillippo example.R

In htx-r/GenericModelNMA:

# libraries
if(!require(dplyr)) {install.packages("dplyr"); library(dplyr)}
if(!require(tidyr)) {install.packages("tidyr"); library(tidyr)}
if(!require(wakefield)) {install.packages("wakefield"); library(wakefield)}
if(!require(ggplot2)) {install.packages("ggplot2"); library(ggplot2)}
if(!require(sandwich)) {install.packages("sandwich"); library(sandwich)}

###########################
#  Simulate AB and AC data
###########################

set.seed(61374988)
# Study characteristics
N_AB <- 500
N_AC <- 300
agerange_AB <- 45:75
agerange_AC <- 45:55
femalepc_AB <- 0.64
femalepc_AC <- 0.8
# Outcome model
b_0 <- 0.85
b_gender <- 0.12
b_age <- 0.05
b_age_trt <- -0.08
b_trt_B <- -2.1
b_trt_C <- -2.5

# AB - IPD

AB.IPD <-
  rbind(
    # Generate A arm
    r_data_frame(n = N_AB/2, # Number of individuals in arm A
                 id, # Unique ID
                 age = age(x = agerange_AB), # Generate ages
                 gender = gender(prob = c(1 - femalepc_AB, femalepc_AB)), # Generate genders
                 trt = "A" # Assign treatment A
    ),
    # Generate B arm
    r_data_frame(n = N_AB/2, # Number of individuals in arm B
                 id, # Unique ID
                 age = age(x = agerange_AB), # Generate ages
                 gender = gender(prob = c(1 - femalepc_AB, femalepc_AB)), # Generate genders
                 trt = "B" # Assign treatment B
    )
  ) %>%
  # Generate outcomes using logistic model
  mutate(
    yprob = 1 / (1 + exp(-(
      b_0 + b_gender * (gender == "Male") + b_age * (age - 40) +
        if_else(trt == "B", b_trt_B + b_age_trt * (age - 40), 0)
    ))),
    y = rbinom(N_AB, 1, yprob)
  ) %>%
  select(-yprob) # Drop the yprob column

# AC - aggregate data
AC.IPD <-
  rbind(
# Generate A arm
r_data_frame(n = N_AC/2, # Number of individuals in arm A
             id, # Unique ID
             age = age(x = agerange_AC), # Generate ages
             gender = gender(prob = c(1 - femalepc_AC, femalepc_AC)), # Generate genders
             trt = "A" # Assign treatment A
),
# Generate C arm
r_data_frame(n = N_AC/2, # Number of individuals in arm C
             id, # Unique ID
             age = age(x = agerange_AC), # Generate ages
             gender = gender(prob = c(1 - femalepc_AC, femalepc_AC)), # Generate genders
             trt = "C" # Assign treatment C
)
) %>%
  # Generate outcomes using logistic model
  mutate(
    yprob = 1 / (1 + exp(-(
      b_0 + b_gender * (gender == "Male") + b_age * (age - 40) +
        if_else(trt == "C", b_trt_C + b_age_trt * (age - 40), 0)
    ))),
    y = rbinom(N_AC, 1, yprob)
  ) %>%
  select(-yprob) # Drop the yprob column


AC.AgD <-
  cbind(
    # Trial level stats: mean and sd of age, number and proportion of males
    summarise(AC.IPD, age.mean = mean(age), age.sd = sd(age),
              N.male = sum(gender=="Male"), prop.male = mean(gender=="Male")),
    # Summary outcomes for A arm
    filter(AC.IPD, trt == "A") %>%
      summarise(y.A.sum = sum(y), y.A.bar = mean(y), N.A = n()),
    # Summary outcomes for C arm
    filter(AC.IPD, trt == "C") %>%
      summarise(y.C.sum = sum(y), y.C.bar = mean(y), N.C = n())
  )


###########################
#  MAIC implementation
###########################

#** 1. Using the method of moments, we estimate the weights in the AB trial:
                                      # logit(w_it) = a_0+a1*X^_EM_it X^Em is the effect modifier- here is age

# the objective function to minimize
objfn <- function(a1, X){
  sum(exp(X %*% a1))
}
gradfn <- function(a1, X){
  colSums(sweep(X, 1, exp(X %*% a1), "*"))
}

# centralize the X's:
# this is done - subtraction: age-age.mean and age^2-(age.mean^2 + age.sd^2) on each row
X.EM.0 <- sweep(with(AB.IPD, cbind(age, age^2)), 2,
                with(AC.AgD, c(age.mean, age.mean^2 + age.sd^2)), '-')
#
print(opt1 <- optim(par = c(0,0), fn = objfn, gr = gradfn, X = X.EM.0, method = "BFGS"))

# ... result
a1 <- opt1$par # the optimal a1
wt <- exp(X.EM.0 %*% a1)

# resclaed weight
wt.rs <- (wt / sum(wt)) * N_AB

summary(wt.rs)
qplot(wt.rs, geom="histogram",
      xlab = "Rescaled weight (multiple of original unit weight)",
      binwidth=0.25)
# the effective sample size is estimated by ..
sum(wt)^2/sum(wt^2)


# the age is balanced now in AB population with AC population after weighting (in terms of mean and SD)
AB.IPD%>%
  mutate(wt) %>%summarise(age.mean = weighted.mean(age, wt),age.sd = sqrt(sum(wt / sum(wt) * (age - age.mean)^2))
)
AC.AgD[,c("age.mean" ,"age.sd")]


#*** 2. estimate the indirect BC effect
# AB effect
fit1 <-
  AB.IPD %>% mutate(y0 = 1 - y, wt = wt) %>%
  glm(cbind(y,y0) ~ trt, data = ., family = binomial, weights = wt)
# Sandwich estimator of variance matrix
V.sw <- vcovHC(fit1)
# The log OR of B vs. A is just the trtB parameter estimate,
# since effect modifiers were centred
print(d.AB.MAIC <- coef(fit1)["trtB"])

print(var.d.AB.MAIC <- V.sw["trtB","trtB"])

# AC parameter
# Estimated log OR of C vs. A from the AC trial
d.AC <- with(AC.AgD, log(y.C.sum * (N.A - y.A.sum) / (y.A.sum * (N.C - y.C.sum))))
var.d.AC <- with(AC.AgD, 1/y.A.sum + 1/(N.A - y.A.sum) + 1/y.C.sum + 1/(N.C - y.C.sum))

# Indirect comparison: BC
print(d.BC.MAIC <- d.AC - d.AB.MAIC)
#
print(var.d.BC.MAIC <- var.d.AC + var.d.AB.MAIC)

###########################
#  STC implementation
###########################

AB.IPD$y0 <- 1 - AB.IPD$y # Add in dummy non-event column
# Fit binomial GLM
STC.GLM <- glm(cbind(y,y0) ~ trt*I(age - AC.AgD$age.mean),
               data = AB.IPD, family = binomial)
summary(STC.GLM)


# Try adding prognostic variables to improve model fit
add1(STC.GLM, ~.+gender, test="Chisq")

# does not improve the fit match so will leave gender out of the model

# estimate the AB effect in the AC population
print(d.AB.STC <- coef(STC.GLM)["trtB"])

# its variance
print(var.d.AB.STC <- vcov(STC.GLM)["trtB","trtB"])

## indirect comparison
print(d.BC.STC <- d.AC - d.AB.STC)
print(var.d.BC.STC <- var.d.AC + var.d.AB.STC)
##################
# Summary
##################

# true AB effect in the AC population: log OR B vs A
d.AB.TRUE <- b_trt_B + b_age_trt * (AC.AgD$age.mean - 40)


# unadjusted estimate of AB effect in AC population
AB.IPD %>% group_by(trt) %>%
  summarise(y.sum = sum(y)) %>%
  spread(trt, y.sum) %>%
  with({
    d.AB.AB <<- log(B * (N_AB/2 - A) / (A * (N_AB/2 - B)))
    var.d.AB.AB <<- 1/B + 1/(N_AB/2 - A) + 1/A + 1/(N_AB/2 - B)
  })

# True BC effect in AC population
d.BC.TRUE <- b_trt_C - b_trt_B

# unadjusted BC estimated effect
d.BC.NAIVE <- d.AC - d.AB.AB
var.d.BC.NAIVE <- var.d.AC + var.d.AB.AB

#
plotdat <- data_frame(
  id = 1:10,
  Comparison = factor(c(rep(1,4), 2, 2, rep(3,4)),
                      labels = c("B vs. A", "C vs. A", "C vs. B")),
  Estimate = c(d.AB.TRUE, d.AB.MAIC, d.AB.STC, d.AB.AB,
               b_trt_C + b_age_trt * (AC.AgD$age.mean - 40), d.AC,
               d.BC.TRUE, d.BC.MAIC, d.BC.STC, d.BC.NAIVE),
  var = c(NA, var.d.AB.MAIC, var.d.AB.STC, var.d.AB.AB,
          NA, var.d.AC,
          NA, var.d.BC.MAIC, var.d.BC.STC, var.d.BC.NAIVE),
  lo = Estimate + qnorm(0.025) * sqrt(var),
  hi = Estimate + qnorm(0.975) * sqrt(var),
  type = c("True", "MAIC", "STC", "Unadjusted",
           "True","Unadjusted",
           "True", "MAIC", "STC", "Unadjusted")
)
ggplot(aes(x = Estimate, y = id, col = type, shape = type), data = plotdat) +
  geom_vline(xintercept = 0, lty = 2) +
  geom_point(size = 2) +
  geom_segment(aes(y = id, yend = id, x = lo, xend = hi), na.rm = TRUE) +
  xlab("Estimate (Log OR)") +
  facet_grid(Comparison~., switch = "y", scales = "free_y", space = "free_y") +
  scale_y_reverse(name = "Comparison in AC population", breaks = NULL, expand = c(0, 0.6))